
import numpy as np
import matplotlib.pyplot as plt

size = 10
nsand = np.zeros([size, size], dtype='uint8')
for i in range(size):
    for j in range(size):
        nsand[i,j] = np.random.randint(0,4)

def avalanch( nsand, size ):
    for i in range( size ):
        for j in range( size ):
            if nsand[i,j] >=4 :
                nsand[i,j] -= 4;
                if i-1 >= 0: nsand[i-1,j] += 1
                if i+1 < size: nsand[i+1,j] += 1
                if j-1 >= 0: nsand[i,j-1] += 1
                if j+1 < size: nsand[i,j+1] += 1
                return 1
    return 0

nstep = 10000
magnitudes = []
for i in range(nstep):
    print("i=",i)
    nsand[(int)(size/2),(int)(size/2)] += 1; magnitude=0
    while avalanch( nsand, size ) > 0 :
        magnitude +=1
        #print("--------------------------------- \n", nsand)
    if magnitude > 0:
        #print("magnitude = ", magnitude)
        magnitudes.append( magnitude )
#print("magnitudes = ", magnitudes)
x = np.bincount(magnitudes)
data = []
for i in range( x.size ): # 去除次数为０的点
    if x[i] > 0: data.append( [ i, x[i] ] )
xdata = []; ydata= [] #
for t in data:
    xdata.append( np.log10(t[0]) ); ydata.append( np.log10(t[1]) )
plt.plot(xdata, ydata, "o")
plt.xlabel(r"$log_{10}$ magnitude")
plt.ylabel(r"$log_{10}$ (number of avalanches)")
plt.savefig("PowerLaw.png")
plt.show()